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SECTION  1 
INTRODUCTION 


The  IvlDCOM  code  is  a digital  computer  program  that  calculates 
ll.l;,  VI. T,  and  IT  electric  and  magnetic  Tield  strengths  in  ambient  and 
nuclear-disturbed  environments.  The  code  is  intended  for  use  when  a 
relatively  detailed  analysis  of  the  propagation  between  two  terminals 
is  required  in  nuclear-disturbed  environments.  The  code  can  be  used 
alone  to  study  the  effects  of  weapon,  atmosphere,  and  propagation 
parameters  on  received  signals  or  can  be  used  in  conjunction  with 
receiver  antenna  and  signal  processing  models  to  evaluate  system  per- 
formance . 

Operational  information  for  UT.DCOM  IV  is  given  in  Volume  1 
(User's  Manual).  A description  of  computer  routines  used  to  model  hl.T 
VI. I , and  IT  propagation  is  given  in  Volume  2 (Computer  Routines). 

This  volume  (Volume  3)  provides  mathematical  descriptions  of  the  pro- 
pagation models  in  a manner  that  does  not  presume  knowledge  of  the 
internal  organization  of  the  UITICOM  code  or  fortran  programming  in 
genera  1 . 

The  ambient  and  nuclear  environment  models  used  in  KhDCOM  are 
taken  directly  from  the  WI.I’H  V 1 code.  These  models  are  described  in 
Reference  1-1  and  are  not  further  discussed  here.  The  propagation 
computational  models  used  in  IVITICOM  IV  include  models  developed  for 
W’l.DCOM  Ml,  models  developed  for  DCOM  (a  code  developed  for  the  MI . I (IN 
office),  and  models  developed  by  the  Naval  Ocean  Systems  Center  (NOSC) 
Descriptions  of  models  developed  for  the  DCOM  code  have  been  taken 
from  the  DCOM  code  documentation  (Reference  1-2).  Description  of 
models  developed  by  NOSC  have  been  obtained  from  References  1-3  and 
1-4. 


Preceding  page  blank 


A description  of  reflection  coefficient  models  used  for  1.1,1  , 

VI.I  , and  Id-  propagation  is  given  in  Section  Section  .3  describes 
the  waveguide  model  used  for  111  and  II  propagation  and  Section  I 
describes  the  skywavc  model  used  for  l.l  propagation.  Introductory 
information  concerning  propagation  model  limitations,  coordinate  sys- 
tems, and  matrix  notation  is  given  below. 

MODEL  LIMITATIONS 

It  should  he  noted  that  there  are  potentially  significant  phys- 
ical processes  that  either  are  not  treated  or  are  very  crudely  modeled, 
since  they  either  have  been  regarded  as  being  of  secondary  importance 
in  past  applications  or  have  not  yet  been  modeled  adequately. 

The  calculation  of  II. I propagation  characteristics  is  difficult 
due  to  the  long  wavelengths  involved.  In  particular,  ionospheric  vari- 
ations transverse  to  the  propagation  path  are  not  modeled  and,  as  a 
result,  predicted  degradation  is  probably  pessimistic.  Standard  IV KB 
approximations,  which  assume  a slow  variation  in  ionospheric  or  ground 
electrical  properties  in  the  direction  of  propagation,  are  used.  In 
addition,  there  is  some  uncertainty  in  selecting  the  appropriate  com- 
binations of  l>- , L-,  and  lower  I -region  ionospheric  models  for  I l.l 
pred i ct i ons . 

Propagation  variations  transverse  to  the  propagation  paths 
are  also  ignored  for  VI. T , but  this  limitation  is  not  as  severe  as  for 
I l.l  because  of  the  relatively  short  VI, I-  wavelengths. 

The  If  model  is  an  approximate  ray  theory  model  using  correc- 
tion factors  for  diffraction  around  the  earth's  bulge.  The  effects 
of  the  geomagnetic  field  are  included  in  the  reflection  coefficient 
calculation,  and  a geometric  search  procedure  accommodates  nonuniformi- 
ties in  the  ionosphere  in  the  direction  of  propagation.  The  model  has 
produced  good  agreement  with  calculations  performed  with  state-of-the- 
art  models  and  assuming  a uniform  ionosphere  for  frequencies  between 
20  and  60  kllz. 
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COORDINATE  SYSTEM  AND  ANTENNA  REFERENCE 


Many  calculations  performed  by  the  system  models  are  performed 
uitli  g comet  ry  references  to  a specified  location,  in  this  report  refer- 
red to  as  a Computation  Reference  Point  (CRPI.  The  Cartesian  coordi- 
nate system  is  the  most  frequently  used  geometry  system.  In  this 
report  , 

x = direction  of  propagation  path  (km) 

>'  7 perpendicular  to  direction  of  propagation  path  (km) 
alt  i tilde  ( km  1 . 

The  origin  is  at  the  earth's  surface. 

The  propagation  models  assume  that  the  radiated  power  is  de- 
fined. for  \ 1 1 and  Id  , total  power  radiated  from  a dipole  with  a 
particular  orientation  may  be  specified.  The  antenna  zenith  angle 
and  orientation  relative  to  the  propagation  path  (p)  are  shown  in 
figure  1-1.  for  I l.l  , the  total  power  and  p arc  specified,  and  y is 
always  i>()  degrees. 
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DIRECTION  OF  PROPAGATION 


Figure  1-1.  Coordinate  system  and  antenna  orientation. 
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SECTION  2 

IONOSPHERIC  REFLECTION  COEFFICIENTS 


INTRODUCTION 

Both  waveguide  and  ray  theory  solutions  require  that  ionospher 
ic  reflection  coefficients  (ratio  of  the  reflected  field  to  the  inci- 
dent field  1 he  computed.  The  ionospheric  reflect  ion  coefficient  ; can 
he  represented  as  a four-component  matrix,  i.e. 


R = 


m«h  i«n 


lit  r'i 


(2-1 


where 

R is  the  reflection  coefficient  matrix 
H means  polarisation  parallel  to  the  plane  of  incidence 
l means  polarization  perpendicular  to  the  plane  of 
i nc  i deuce . 

The  first  subscript  applies  to  the  polarization  of  the  incident  wave, 
and  the  second  subscript  applies  to  the  polarization  of  the  reflected 
wave . 

The  couplin'.1  coefficients,  R and  R , niav  he  as  large  or 
1 x ii  x 

larger  than  the  primarv  coefficients  R and  R for  normal  nighttime, 
but  they  are  typically  much  smaller  than  the  primary  coefficients  for 
daytime  conditions,  lor  a normal  ionosphere,  the  coefficients  are  a 
function  of  propagation  azimuth  and  the  earth's  magnetic  field.  A 
nuclear  disturbance  results  in  lowering  the  reflection  region;  as  a 
consequence  of  the  increased  electron-neutral  collision  frequency  in 
the  reflecting  region,  a strong  disturbance  can  cause  the  coupling 
coefficients  to  be  essential lv  zero. 
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II  and  VIT  waves  are  reflected  from  the  l>  region.  Most  of 
the  reflected  energy  is  returned  from  the  50-70  km  altitude  level  for 
normal  daytime  ionospheres  and  from  80-00  km  altitude  level  for  normal 
ni  *ht t i nu  ionospheres.  In  these  altitude  ranges,  particularly  for 
ni  ;httime  conditions,  the  electron-neutral  collision  t requcncy  is  com- 
parable to  the  electron  gyrof rcquency , and  the  earth's  magnetic  field 
has  a significant  effect  on  the  ret  lection  coefl icient. 

11,1  reflection  coefficients  are  affected  by  ionisation  in  the 
I and  I -regions  in  addition  to  ionisation  in  the  0- region.  I'htis , 
effect:  from  the  earth's  magnetic  field  are  always  important.  I he  ion 
concentrations  from  ground  level  up  to  the  D-region  also  affect  the 
1,11  reflection  coefficient  values. 

The  hi  DCOM  code  employs  two  calculation  procedures  for  deter- 
mining the  ionospheric  reflection  coefficients.  fhe  first  method  is 
an  iterative  integration  of  a first-order  differential  equation  (Refer- 
ence 7-1).  It  is  applicable  to  both  isotropic  and  anisotropic  iono- 
spheres. The  second  procedure  is  an  application  of  a simple  recursive 
technique  (Reference  7-2)  and  is  used  only  with  isotropic  ionospheres. 

Both  reflection  coefficient  computation  procedures  allow  for 
arbitrary  vertical  profiles  for  electron  and  ion  densities.  These 
profiles  are  divided  into  horizontal  slabs.  The  iterative  integration 
method  allows  a simple  variation  of  the  the  slab  electromagnetic  pro- 
perties. This  procedure  is  used  for  all  ITT  cases  and  tor  Vl.F  and  IT 
cases  where  the  ionosphere  is  highly  disturbed.  I he  recursive  technique 
requires  less  computation  time  than  the  iterative  integration  method. 

COMPUTATION  OF  REFLECTION  COEFFICIENTS 
FOR  ANISOTROPIC  IONOSPHERE 

the  reflection  coefficient  calculation  technique  described  by 
Budden  (Reference  2-1)  is  used  here.  The  model  is  patterned  after  that 
described  by  Sheddy  et  al  (References  2-.S  and  2-4). 

All  ionospheric  and  geomagnetic-dependent  inputs  are  calculated 
prior  to  execution  of  the  model.  A fictit  ions  linear  variation  of  the 
permittivity  with  altitude  is  used  to  simulate  the  earth  curvature. 
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Basic  Equations  Defining  R 

The  ionospheric  reflection  coefficient  matrix  defined  in  liqua- 
tion 2-1  is  found  by  the  numerical  integration  of  differential  equations 
o f t he  ma  t r i x f o nil 


2 1 ll'  - S_m  + S,0R  - R S_n  - R S,  ,R  , ( 2-2 j 

where  the  prime  represents  the  derivative  with  respect  to  the  normal- 
ized parameter  kz  , and 

k = wavenumber  (km  | 
z = altitude  (km) 

The  elements  of  matrix  S arc  partitioned  for  convenience,  with 
the  partitioned  values  being 


( S 1 la  + S1  lb)  ~S1. 


(-dlla+dllb)  ~S 1 2 


(2 -3a) 


( 2 -5b  I 


( d 1 1 a - dl  lb)  dl 


5 ) > 


IS..  - S . J d,. 


I la  1 I hi 


J-,, 


-S. 


(2-3d) 


The  elements  of  the  partitioned  S matrix  are  functions  of  the 
ionospheric  properties,  the  geomagnetic  field  strength,  the  wave  fre- 
quency, and  the  angle  of  incidence  of  the  wave  with  the  ionosphere. 
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Those  matrix  elements  are  defined  in  Appendix  A in  terms  of  an  inter- 
mediate matrix  T and  the  susceptibility  matrix  M. 

Considerable  simplification  occurs  in  the  S,  I,  and  M matrices 
when  the  geomagnetic  field  has  a negligible  influence  on  the  suscepti- 
bility matrix.  In  the  KflK'OM  code  it  is  assumed  that  the  altitude 
range  where  the  magnetic  field  influence  is  negligible  where 


A 


• 111 


(2-1) 


. = tree-space  permeability  (henry/m) 

e = electron  charge  (coulomb) 

II  static  magnetic  field  strength  (oersted) 
m = electron  mass  ( kgj 

" frequency  (radians  s ') 

= electron-neutral  collision  frequency  (s  *) 

j = . 

lhe  lower  altitude  is  determined  entirely  by  electrons.  The  contri- 
bution of  V of  the  normal  positive  and  negative  ion  concentrations  is 
much  smaller  than  the  electron  component. 

VI. 1-  AND  II  OASIS.  Ixpanding  liquation  2-2,  equations  are  writ- 
ten for  the  four  components  of  the  reflection  coefficients  as  they  arc 
computed  for  VI. I and  Id  cases: 


! 


r — - 


->A  - [sna-snh*  At‘l,hraiibl  - s” 

* A J2I  * d22  A]."  * [d12(1*A) 

* A si:"\v] 

->&  = [-K..‘dlli.  -JUh>  * S12  J>  * A J22 

- |Slla  *Sllb’  - A:]  \ * [Ai'1*/, 

+ R d . . ( 1 + R )| 

21  i.  j.  J 

2 1 R'  = R - 2S  + R d_,  + S.  , 111  R 

• l l L 22  jl  l 22  i 21  12  j. J j.  j. 

+ I’d.,,  + R R (d.  . - d,  ) 

122  j.  n ii  l 11a  lib 

+ R S + R d 1 
j.  12  i 2 1 J 


(dlla  ~ d 1 1 b 1 


( 2 - 5b  | 


(2-5d ) 


where  j = \J- 1 

1 LI  CASTS.  Numerical  difficulties  can  occur  when  using  liqua- 
tion 2-5  for  r; LI  reflection  coefficient  calculation.  The  numerical 
difficulties  result  from  short  computer  word  length.  Using  double 
precision  techniques  or  computers  possessing  larger  word  lengths  elim- 
inates the  problem.  A separate  algorithm  was  formulated  for  Kf DCOM 
for  computational  efficiency  and  adaptability  to  machine  type.  The 
scaling  and  factorization  of  the  reflection  coefficient  equations  arc 
based  on  the  following  observations: 

1.  The  S matrix  elements  are  large  (of  order  10(>)  at  I . I . I . 


S 1 1 b J *dllb  'J  U41 


S12  ^ d 1 2 


S21  ^ 'd21  ^ 1 5 1 


s„  **  ~ T,,/C 


where  (!  is  the  cosine  of  the  angle  of  incidence. 

2.  The  diagonal  elements  of  the  R matrix  are  character- 
i st ica 1 1 v 


R ~ - R ~ 1 , (2-6) 

II  ll  _L  J_ 

and  the  coupling  elements  R and  R arc  much  smaller 

X il  i x 

than  unity. 

Ihese  observations  show  that  subtraction  of  large,  nearly  equal  values 
occurs  frequently  in  liquation  2-5.  They  also  provide  the  basis  for 
modification  of  liquation  2-5  for  computational  efficiency. 

lor  1.1,1  cases,  liquation  2-5  is  modified  so  that  the  derivative 
of  the  logarithm  of  R is  integrated,  where 

f{n  R)'  = | (2-7) 

K 

liquation  2-5  is  further  modified  employing  the  observation  noted  in 
Iquation  2-6.  first,  the  following  relationships  arc  made: 

R = 1 - A,  (2-8) 

ii  ii  ] ”i 


xKx  = A2  - 1 • (2-9) 

Substitution  of  liquations  2-8  and  2-9  is  made  into  liquation  2-5.  After 
simplification  of  terms. 
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The  "T"  elements  in  liquations  2-10  belong  to  an  intermediate  matrix 
which  relates  the  susceptibility  matrix,  M,  to  the  S matrix.  The  T 
and  M matrices  are  defined  in  Appendix  A.  Implementation  of  liquation 
2-10  for  iil.r  cases  removes  the  numerical  difficulties  described  earlier 

Integration  Technique 

The  integration  technique  is  that  defined  in  Reference  2-3. 

It  is  performed  using  a fourth-order  Runge-kutta  technique.  The  inte- 
grand after  n+1  steps  is  defined  by 
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Integration  of  liquat  ion  2-5  or  2-10  is  carried  out  In  starting 

at  a high  enough  altitude  (defined  as  ! and  integrating  downward, 

each  step  employing  the  control  criteria  defined  above.  Integration 

is  terminated  based  on  a criterion  which  indicates  that  a free -space 

assumption  can  apply  (defined  as  2,  ). 

hot 

During  the  integration,  the  step  size  is  controlled  by  the 
procedure  described  above!  However,  the  step  size  is  never  allowed 
to  be  less  than  A: 


111 1 n 


A suitable  value  of  \z  . was  found  to  be: 

in  1 n 


in  1 n 
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where  is  the  frco-space  wavelength. 

Initial  Values  for  R 

The  altitude  z^  representing  the  ionosphere  top  is  taken 

high  enough  so  that  reflection  front  higher  altitudes  is  negligible. 

The  ionosphere  is  approximated  as  a uniform  semi-infinite  anisotropic 

slab  above  the  top  altitude,  and  thus  only  upgoing  waves  exist.  The 

initial  step  is  to  determine  R at  z . The  procedure  used  here  is 

top  1 

that  of  Shoddy  (Reference  2-5),  and  the  details  are  included  in 
■\ppcndix  B. 


Initial  Values  of  z.  and  z, 

top  bot 

The  integration  starting  altitude  z has  a strong  influence 

s top  c 

on  the  computational  efficiency.  Factors  of  2 to  10  in  computational 
efficiency  can  be  gained,  without  sacrificing  accuracy,  by  the  selec- 
tion of  the  proper  starting  altitude. 

The  electromagnetic  wave  behavior  in  free-space  can  be  treated 
analytically.  Thus,  additional  computational  efficiency  is  gained  by 
defining  the  integration  stopping  altitude,  Z|  ^ , as  the  highest  alti- 
tude where  the  free-space  assumption  can  be  made. 


VI. F AND  IT  OASIS.  The  hi  DCOM  D- region  ionization  profile  is 
defined  up  to  100  km.  A satisfactory  starting  altitude  criterion  ap- 
plicable for  all  undisturbed  profiles  could  not  be  determined.  However, 
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The  ionosphere  top  for  VI.I-/I.I  eases  is  defined  as  that  a It  i 
tude  above  which 


B ' 2 maximum 
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The  ionosphere  bottom  is  defined  as  that  altitude  below  whieh 


B ■ 0.0002 2 . 
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l id  ll.ASliS.  Ions  at  altitudes  are  important  at  I 1. 1 , and  the 

re  bottom  altitude,  , is  s 

bot 

t ion  of  the  upper  (start i ng ) altitude. 


ionosphere  bottom  altitude,  z,  , is  set  to  '.’round  level.  Determina- 

bot 


top 


is  more  complex,  be 1 ini ng 


t op 


too  high  increases  the  computation  time  considerably  without  add- 


accuracy.  Moreover,  the  algorithm  used  in  YI.F/IJ  cases  for  deter- 


mining .■  has  little  value  here,  as  it  is  more  applicable  for  D-  rev  i on 

t op  11 


effects  and  not  for  the  I - and  1 -regions  where  magnetic  field  effects 
are  strong.  Hie  algorithm  used  in  IVITKIOM  uses  an  anisotropic  defini- 
tion of  the  squared  index  of  refraction,  liquations  2-l(>,  where  for  the 
ordinary  wave  (see  Reference  2-<>) 
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The  UT.OCOM  1.1.1  model  defines 
sion  loss  due  to  nondeviate 


: as  the  altitude  where  the  transmi 

t op 

incremental  absorption  satisfies 
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(2-19) 


i = 8.(>8h  ' Imln) 

I c 

c = velocity  of  light  (kin/sj 
Im(n)  - imaginary  part  of  the  index  of  refraction. 


If  liquation  2-19  is  not  satisfied  up  to  200  km,  then 
the  value  of  200  km. 


top 


is  given 


REFLECTION  COEFFICIENTS  FOR 
ISOTROPIC  IONOSPHERE 


General 

The  reflection  coefficients  for  an  isotropic  ionosphere  arc 
calculated  using  a method  described  by  Wait  (Reference  2-2).  liarth 
curvature  and  magnetic  field  effects  are  ignored.  The  cross  terms  of 
R I R and  R ) are  zero,  and  the  diagonal  elements  ( R and  R ) are 
calculated  independently  of  one  another.  The  following  relationships 
are  defined  for  convenience  to  help  avoid  duplication  of  similar 
equat ions : 
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The  ionized  layer  is  represented  by  homogeneous  slabs  as  shown 
n figure  2-1,  where  the  wave  is  incident  on  the  layer  bottom  by  angle 
. and 


Figure  2-1.  Plane  wave  incident  on  stratified  medium 
with  homogeneous  slabs. 


n = refractive  index  of  mth  slab 
in 

n = refractive  index  of  free  space  = unity 
o 

Az  = thickness  of  mth  slab 
m 

z = altitude  above  which  there  arc  no  downcoining  waves 
top 

z.  = altitude  at  which  profile  is  truncated, 
hot 


I'he  region  below  the  stratified  medium  is  considered  to  be  free  space. 


B.v  imposing  the  condition  that  the  Mth  slab  is  semi  finite, 
onl>  upgoing  waves  are  allowable  in  the  Mth  slab  and  a recursive  solu- 
tion for  the  reflection  coefficient  is  obtained. 

lor  either  a vertically  or  horizontally  polarized  plane  wave 
incident  on  the  bottom  ot  a stratified  medium,  the  reflection  coeffi- 
cient (ratio  of  downgoing  to  upgoing  wave  amplitude)  is  defined  as 


v,h  K + I 
o ] 


(2-21  ) 


where  is  the  zeroth  term  in  a series  of  K^'s  and  Z is  the  first 

term  in  a series  of  Z 's. 

in 

1 s ;l  function  ot  the  squared  index  of  refraction  of  the 

mth  slab: 
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tor  vertical  polarization 
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for  horizontal  polarization 


b = \ /s  in“  ' - ’ ~ 

m c \/  i m 


(2-23) 


KC)>  :l  characteristic  impedance  of  free  space,  is  simply  defined  as 


COS  P . 


~ — for  vertical  polarization 
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~~  for  horizontal  polarization 
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(2-2  5d  ) 


in  the  above  equations, 

f = permittivity  of  free  space  (farads  m ') 
c = velocity  of  light  (m  s') 

The  solution  is  obtained  by  starting  at  the  Mth  slab,  computing 
and  K I , from  them  computing  Z,.  j (using  liquation  2-25d  to  define 
Z ),  and  then  proceeding  downward,  repeating  the  procedure  (using  liqua- 
tions 2-25a,  2-25b,  and  2-25c)  until  a suitable  point  for  truncating 
t he  stratified  medium  is  reached. 

Analytic  Representation  of  Isotropic 
Reflection  Coefficient 

Depressed  ionospheres  produce  reflection  coefficients  with 
little  magnetic  field  dependence.  Moreover,  earth  curvature  effects 
are  reduced.  An  analytic  representation  of  the  isotropic  reflection 
coefficient  is  applicable  for  these  cases,  and  its  use  can  result  in 
considerable  saving  in  computation  time. 

This  approximate  formulation  is  a modified  form  of  that  sug- 
gested by  Wait  (References  2-7  and  2-8).  It  is  especially  convenient 
for  use  in  the  VLF  mode  solutions,  which  require  an  analytical  expres- 
sion for  the  reflection  coefficient. 

The  reflection  coefficients  arc  calculated  for  fixed  angles 
of  incidence  and  are  represented  by 
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whore  ('  is  the  cosine  of  the  angle  of  Incidence  at  the  ionosphere  and 
the  's  are  determined  as  follows.  Reflection  coefficients  are  calcu- 
lated for  two  real  angles  of  incidence  with  cosines  (i*  and  . If 


!•.((:.)-=  jrj  exp  ( j 


r . ( C -, ) = |r^|  exp 


where  and  I are  the  phases  (radians)  referred  to  the  same  reference 
a 1 1 i tude  , z , t hen 
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The  phase  of  the  isotropic  reflection  coefficient  is  referenced 
to  the  bottom  of  the  ionosphere.  To  use  the  analytical  reflection  coef- 
ficient representation  just  described,  the  phase  must  he  referenced  to 
the  same  altitude  for  all  incident  angles  of  interest.  The  phase  refer- 
ence altitude  is  determined  from  the  altitude-dependent  refractive  index. 


r.(C  ) = -exp  (itjl.  ) exp 
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(2-26) 
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When  the  altitude  variation  of  n“(z)  is  dominated  by  variations 
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in  B ( z ) (see  liquation  2-10),  reflection  maximizes  (Reference  2-‘>|  at 
approximately  the  altitude  z , where 

B(zr)  = 2 cos2  4>.  , (2-34) 

where  <j>.  is  the  angle  of  incidence  at  t he  ionosphere. 

A fixed  reference  altitude  is  required,  and  a suitable  value 
is  defined  by 


B(zr)  = 0.04 


by  noting  that  cos  <Ji  for  angles  of  interest  is  near  the  range  0.1  to 

0.2  . 

The  exact  altitude  used  for  the  phase  reference  altitude  is 
not  critical.  The  calculation  procedures  are  defined  so  that  varia- 
tions in  ZR  that  arc  not  large  relative  to  z have  little  effect  on 
the  final  answer.  The  phase  of  the  reflection  coefficient,  referred 
to  the  reference  altitude,  is 
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(2-36) 


where 


)>R  = phase  referred  to  reference  altitude  (radians) 

= altitude  where  the  recursive  solution  is  stopped 

(km) 

(f)  . = phase  at  z.  . returned  from  recursive  calculations 

m i n 1 hot 

(rad i ans ) . 
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SECTION  3 

WAVEGUIDE  SOLUTIONS 


INTRODUCTION 

The  LLL/YLI  propagation  models  compute  the  received  electric 
field  strength  ;is  the  summation  of  waveguide  modes.  Mode  solutions 
are  found  for  a uniform  waveguide.  In  a strict  sense,  this  solution 
is  only  applicable  for  I l.l/VT.I  transmissions  between  a uniform  concen- 
tric ionosphere  and  earth  with  a uniform  magnetic  field.  In  situations 
where  the  ionosphere  is  disturbed  by  nuclear  detonations  or  where  the 
propagation  path  crosses  the  sunset  or  sunrise  terminator  or  a land- 
sea  boundary,  the  waveguide  is  not  uniform  along  the  path. 

Two  procedures  are  usid  in  IVLDCOM  to  account  for  this  nonuni- 
formitv.  The  first  method  is  an  approximation  which  assumes  that  the 
ionospheric  or  ground  properties  vary  slowly  in  the  direction  ot  pro- 
pagation. This  method  is  called  the  KkB  approx imat ion.  The  second 
method  models  the  nonuniform  propagation  path  as  several  adjacent  uni- 
form waveguides  with  different  propagation  characteristics.  An  energy 
conversion  process  (mode  conversion)  is  computed  at  each  interface. 

The  Mode  Conversion  method  is  used  only  for  ALT  cases.  The  WKB  method 
is  an  option  available  for  ALT'  cases  and  is  used  for  all  ILL  cases. 

The  computational  sequence  requires  defining  uniform  waveguide 
solutions  (eigenvalues),  using  these  to  define  values  for  attenuation 
rates  and  height  gain  and  mode  excitation  factors,  and  applying  these 
values  in  mode-summing  equations  to  compute  total  field  strength.  The 
procedures  arc  described  in  order  below. 

THE  GENERAL  MODE  EQUATION 

The  general  mode  equation  for  a uniform  waveguide  can  be  writ- 
ten in  the  form 
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where  R , R , R , and  R are  the  ionosphere  reflection  coefficients 

II  II  ’ x x X II  II  x 

described  in  Section  2,  and  MR|(  and  R are  ground  reflection  coeffi- 
cients for  vertically  and  horizontally  polarized  waves.  The  mode  equa- 
tion requires  that  the  R's  and  R's  can  be  referenced  to  the  same  alti- 
tude. The  R's  are  referenced  to  the  bottom  of  the  ionosphere,  or  to 
an  altitude  where  most  of  the  energy  is  reflected,  or  to  ground,  de- 
pending on  the  computational  option  used,  liquations  for  the  It's, 
referenced  to  a specified  altitude,  were  taken  from  Reference  3-1  and 
adapted  for  use  with  Airy  functions.  The  equations  are  described  in 
detail  in  Appendix  C. 

VLF  EIGENANGLE  DETERMINATION 
General 

The  e i genanr 1 es , " , are  those  values  for  which  I I1')  vanishes, 
n 

They  are  difficult  to  determine  explicitly  and  are  found  by  iteration 
techniques.  Determination  of  •'  in  Ui  DCOM  is  accomplished  by  two  pro- 
cedures. The  first  (identified  here  by  A)  is  a search  technique  over 
the  complex  ’-plane,  and  is  an  adaptation  of  the  MODI  SROII  algorithm 
i Reference  3-1  I.  The  second  procedure  (identified  here  by  Bl  employs 
an  analytic  approximation  of  F(0). 

Both  procedures  use  an  analytic  representation  of  the  ionospheric 
reflection  coefficients.  This  can  be  illustrated  with  the  aid  of  Figure 
3-1.  Procedure  \ first  identifies  a set  of  rectangular  regions  of  the 
-plane  where  the  anisotropic  ionospheric  reflection  coefficients  are 
computed  at  each  corner.  A l.agrangian  interpolation  representation 
of  the  reflection  coefficients  is  made  for  a rectangle  followed  by  a 
l urch  method  to  locate  the  eigenangles.  For  the  purposes  of  KF.DCOM, 
these  approximate  eigenangles  have  been  found  sufficiently  accurate 
and  iteration  to  the  exact  value  is  not  performed  as  in  the  M0DRSRC1I 
program. 

Procedure  B computes  the  isotropic  ionospheric  reflection  co- 
efficients at  two  real  0 values.  An  analytic  formulation  of  the 
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Figure  3-1.  The  complex  0-plane  illustrating  the  6 values  used  in 
the  ionospheric  reflection  coefficient  computations 
for  initial  eigenangle  determination. 

reflection  coefficients  (using  the  tw.o  precomputed  values)  is  employed 
to  estimate  approximate  eigenangles.  Next,  a Ncwton-Raphson  algorithm 
with  these  approximate  eigenangles  as  initial  values  iterates  to  an 
exact  eigenangle  value  using  the  anisotropic  reflection  coefficient 
model.  This  procedure  is  described  in  detail  later  in  this  section. 

Procedure  B is  fast  and  accurate  for  daytime  or  strong  nuclear 
depressed  ionospheres.  The  approximate  eigenangles  are  sufficiently 
accurate  for  severely  depressed  ionospheres,  and  the  exact  eigenangle 
determination  is  not  made.  Procedure  A utilizes  more  computation 
time,  hut  is  much  more  accurate  for  weakly  disturbed  nighttime  iono- 
spheric profiles.  To  help  illustrate  the  weakness  of  Procedure  B, 
Figure  3-2  shows  the  loci  of  eigenangles  as  a daytime  profile  is 
gradually  changed  to  a nighttime  profile.  The  eigenangles  start  at 
the  left  for  the  Tl  and  I'M  modes,  and  progress  as  indicated  by  the 
arrow  as  the  ionosphere  is  varied.  For  Procedure  B,  note  that  the 


Loci  of  eigenangles. 


isotropic  solution  provides  an  adequate  approximation  for  the  more 
depressed  ionospheres,  and  convergence  to  the  exact  modes  will  occur 
with  some  confidence.  However,  as  the  nighttime  profile  is  approached 
the  eigenangles  for  the  first  IT.  and  'I'M  modes  arc  nearly  equal,  and 
Procedure  B will  result  in  wrong  convergence.  The  Search  Method  (Pro- 
cedure A)  works  well  for  the  nighttime  cases,  even  those  where  the 
eigenangles  are  almost  identical.  Both  procedures  will  now  he  des- 
scribed  in  more  detail. 

Search  Procedure  A 

\ s mentioned  above,  this  procedure  is  an  adaptation  of  the 
MOPbSRCII  techniques.  A complete  description  is  found  in  Reference  5-1, 
hut  the  procedure  is  outlined  here  to  assist  in  the  description  of 
the  U'hDCOM  models  (Volume  2).  The  search  procedure  is  identical  to 
all  rectangles,  l-.ach  rectangle  is  5 degrees  wide  along  the  Re(0)  axis. 
The  length  along  the  Imp)]  axis  is  frequency  dependent  and  is  computed 
from  a preselected  maximum  attenuation  rate  as  a function  of  Re ( 0 ) 

(see  figure  5-1). 

Die  eigcnangle  " search  for  verv  oblique  ionospheric  incidence 
1 ‘ n • 

is  numerically  difficult  due  to  a pole  in  l:(b)  at  0 = 90°,  and  another 
potential  pole  in  the  ground  reflection  coefficients,  Rn  and  R . 

Thus,  the  modal  equation,  1 •'  ( 1 ' ) , is  modified  by  first  defining 
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where'  F is  the  cosine  of  the  ancle  of  incidence,  and  N , , I > , \ , 

and  I)  are  defined  with  no  poles  in  the  region  of  interest.  Ihe 
pole  at  ■ = 90°  has  been  removed  from  the  U's  with  the  definition  of 
the  X's  (Reference  5-2).  A new  modal  equation  is  formulated  with  the 
substitution  of  liquations  5-2  through  5-7  into  liquation  5-1  and  re- 
arranging: 


I)  D 
II  n x x 


I (9) 


(5-8) 


The  modified  modal  equation  F (n)  is  analvtic  and  presents  little 
numerical  difficulty  for  the  1 region  of  interest. 

The  validity  of  the  search  procedure  is  premised  or  the  Argu- 
ment Principal  of  complex  variable  theory,  which  insures  that  the  phase 
change  of  1 (')  will  complete  a cvcle  for  each  eigenangle  when  I (')) 
is  traced  counterclockwise  around  the  rectangle. 

Figure  5-5  presents  an  enlargement  of  one  of  the  rectangles 
of  Figure  5-1.  The  finer  mesh  site  is  smaller  than  the  separation  of 
the  eigenangles.  Starting  at  the  rectangle  corner  A.  and  proceeding 

counterclockwise,  the  character  of  1 I")  is  examined  at  each  mesh  cor- 

m 

ner.  This  investigation  is  based  on  the  following  definition  of  F (0): 

m 


i . 


= >R  + -i,: 


(5-9) 


is,  when  i (9)  = 0,  1,  I.  ‘0,  and  ' is  undefined.  Also  note 
m R I m 


other  I (9)  characteristics  in  Table  5-1. 
m 


rhus 


i f 


changes  sign 


from  one  mesh  corner  to  an  adjacent  corner,  then  the  line  F = 0 enters 
that  mesh  somewhere  along  the  mesh  edge.  Thus,  using  the  example  il- 
lustrated in  Figure  5-2,  no  F sign  change  would  occur  until  mesh  5. 

The  eigenangle  search  continues  as  follows: 
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1.  The  mesh  3 edge  where  I:  =0  exist  is  determined  by  com- 
puting F ;it  the  remaining  two  corners.  A check  is  made 
to  determine  if  possibly  one  other  F = 0 curve  may  also 
enter  and  exit  mesh  3 (all  adjacent  corners  would  have 
an  opposite  sign).  Finally,  F is  computed  at  the  four 
corners  to  see  if  an  F = 0 line  enters  and  exits  mesh  3. 
These  tests  assume  that  F j and  F^  are  linear  with  0 
along  the  edges,  and  so  only  one  Fj  = 0 or  F = 0 line 
crossing  can  occur  along  one  edge. 

2.  The  F = 0 line  exits  along  the  upper  edge  of  mesh  3. 

Mesh  2 is  examined  as  in  step  1 followed  by  meshes  7,  (>, 

and  11.  The  exit  point  from  the  rectangle  is  noted  (mesh 
11)  so  that  reentry  will  not  be  performed  as  the  search 
continues  around  the  rectangle  contour. 

3.  Fines  of  Fj  0 and  F = 0 were  found  in  mesh  7.  Again, 

assuming  the  linear  behavior  of  F^  and  F along  the  mesh 

edges,  simple  curves  (Hyperbolic!  are  formed  within  the 
mesh.  If  the  computed  crossing  lies  along  the  F = 0 line 
being  traced,  an  eigenangle  is  identified.  The  remaining 
rectangle  is  similarly  searched,  starting  with  mesh  1 and 
continuing  until  all  four  rectangle  edges  have  been  traced 

The  cigenangles  found  arc  approximate,  as  the  mesh  employed  a 
simple  hyperbolic  curve  mesh  geometry  to  determine  the  eigenangle  lo- 
cation. A Newton-Raphson  iteration  procedure  is  used  to  remove  this 
approximation.  When  two  cigenangles  are  located  very  close,  this 
procedure  is  modified  to  insure  that  the  Newton-Raphson  procedure  does 
not  iterate  both  approximations  to  the  same  eigenangle. 

Isotropic  Ionosphere  Procedure  B 

Procedure  II  initially  assumes  an  isotropic  ionosphere,  which 
is  applicable  when  the  earth's  magnetic  field  is  small  and  R and 

ii  f,  x M 

R are  approximately  zero.  For  this  situation  liquation  3-1  can  be 
rewritten  as 
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1(0}  = (1  = R R ) ( 1 - R R ) 

X XX  X 


(3-10} 


ihis  permits  independent  solutions  for  vertical  and  horizontal  polari 
cation  to  he  made;  ie, 


I - R R = 0 


1 = R R =0 

X 11  X 


13-11} 

(3-12) 


lor  convenience  in  the  following  discussion,  the  definitions 
he  low  are  made: 


R 

II  II 


i A 


R 

ii  n 


when  solutions  for  vertical  polari- 
zat ion  arc  sought 

(3-13 i 

when  solutions  for  horizontal 
po 1 a r i zat i on  a re  sought 


when  solutions  for  vertical  polari- 
zation arc  sought 

(3-14) 

when  solutions  for  horizontal 
polarization  arc  sought 


I.quations  3-11  and  3-12  can  he  written  in  terms  of  Airv  func- 


t ions  as 


Alt  ) B(t  ) = e 
n n 


i 2 rrn 


1 ••>  - 1 i>  i 


lor  the  nth  mode  and  where  t is  a uniquely  determined  eigcnangle 
(Reference  3-2).  The  functions  A ( t ) and  B(t  ) are  defined  hv 


A'V  ^ 


Bit 


V'n  - -VR 1 + W*,,  - -V 

W2(tn  = ^R1  + l«iW2(tn  ‘ > R 1 
>V2(tn'  ~ (lU2(tn) 

W!(V  - ('W 


(3-16} 


(3-17) 
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where  U ^ ( t ) and  IV  t ) are  Airy  functions  of  the  first  and  second  kinds 
and  the  prime  indicates  the  derivative  with  respect  to  the  argument. 

In  liquations  3-10  and  3-17,  q.  and  q are  the  normalized  imped- 
ance is  defined  by 


T = 


1 /2 

i (y„  - t ) ' ( 1 - r. 
• ' R n i 

1 + r. 

i 


(3-18) 


where 


H k "R 
a 


{3-19} 


z = reference  altitude  (km),  defined  by  l.quation  2-35 

a = earth's  radius  (km) 

k = wavenumber  (km) 

k = (ka/2) 1/3 
a 

The  normalized  surface  impedance  at  the  ground  is 


jvni 

/ - i k n 
V ’ a g 


vert i ca 1 po 1 ar i cat  ion 


hor i zonta 1 po 1 a r i zat i on 


(3-20) 


where 

n"  = > - i 

U rot. 

o 

2 2 1/2 
Q = j k (n  - .3  ) 

II  ' a g 

S = sine  of  the  cigenangle  at  the  ground 
k = wavenumber  (km  ) 

= relative  dielectric  constant  of  the  earth 
o = earth  conductivity  (mhos/m) 
t = permittivity  of  free  space  (farads/m)  . 

The  above  definitions  coupled  with  the  analytic  approximations 
to  the  reflection  coefficients  (liquation  2-26)  permit  several  approxi- 
mations to  liquation  3-15.  liquation  3-15  can  be  approximated  by  (Refer- 
ence 3-2) 
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1 

j 


whe  re 


- j2  n 
c • 


r.  r I 

■ 8 8 


oxp 


f 3-21  ) 

(3-22) 

(3-23) 


<:  is  the  complex  cosine  of  the  angle  of 
orul  C is  the  complex  cosine  of  the  angle 
Mso,  c and  are  related  In 


incidence  at  the  ionosphere, 
of  incidence  at  the  ground. 


!•  i s a 
8 

of  Hit) 


mode  as 


C = (<!' 


correction  factor  to  account 
and  is  discussed  below.* 

I*s  1 n8  I quat  ion  2-2(' , l.puat  ion 


(3-24 ) 

tor  errors  in  the  approximation 
■■>-21  can  he  written  for  the  nth 


+ ' i *|i:'  + J ‘ ) - <-(<:'  j + fee' ) 

+ J n ( r ) + i mi  1 1-  ) 

8 " 8 


3 25 


If  rg  is  represented  by  a Iresncl  reflection  coefficient, 
expressed  as 


can  be 

8 


h ' ( t ) - (|i\  it) 
u j ( t ) - qW j It  ) 

j exp  r.i 4/5 1 -t  )'v  -1  (JJ—LLL I 

L J ‘ i ( - 1 ) 1 ' - / 


* 


for  simplification  the  subscript  "n"  has  been  omitted  in  this 
ur  * n • If  Mil  he  understood  that  t refers  to  the  nth  mode 
either  vertical 'or  horizontal  polarization. 


sect  Lon 
for 
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For  vortical  polari  at  ion,  1 varies  rapidly  with  ..i.und  conduct  iv  i t> 
i ■’  the  conduct  iv  it\  is  loss  than  10  mho/in. 

AI’PUOA  ! MAT  I ONS  10  ISOTROPIC  MODI;  FQIIAT  I ON  . lo  solve  l.quation 

a 'a,  it  is  convenient  to  define  the  product  I r as 

S g 

h , 1 1 ) qlV,(tl 
i 

iv'lt  ) - qlV  It  l 

1-  , I r - -pp.  . (3-27) 

" K exp^j  1/3  ( -t  l'1  “ 

The  procedure  used  to  solve  l.quation  3-2S  is  as  follows: 


1.  Assume  that  t lies  between  some  limiting  values  for  which 
an  approximate  form  for  the  Airy  function  can  he  used. 

2.  Approximate  F hv  its  value  for  highlv  conduct  ini'  earth, 
and  obtain  a solution  to  I . qua t ion  3-25. 

3.  lest  to  see  if  the  solution  for  t is  in  the  assumed  range. 
If  not,  assume  another  range  and  repeat  Steps  1 and  2. 

1.  It'  the  actual  ground  conductivity  is  less  than  10  mho/m, 
I is  redefined  bv  substituting  the  solution  to  Fount  ion 

;;r 

5-21  from  Step  2 and  using  a value  of  condictivity  slight- 
ly less  than  10  mho/nt.  This  process  is  continued,  re- 
ducing the  conductivity  each  time,  until  the  conductivity 
is  reduced  to  the  desired  low  value. 


No  iteration  on  conductivity  is  required  for  horizontal  polari- 
sation, since  the  high-conductivity  approximation  applies  for  all  rea- 
sonable values  of  ground  conduct i v i tv.  Approximate  values  for  I and 

ii  gr 

the  resulting  solutions  for  liquation  5-25  for  both  vertical  and  hori- 
zontal polarization  are  given  below. 


base  la.  Re(-t)  1,  vertical  polarization. 

This  is  the  modified  flat-earth  approximation,  which  is  better 
when  used  with  waveguide  upper  boundaries  typical  of  normal  daytime 
or  nuclear-disturbed  ionospheres  anil  for  higher  modes,  for  this  case, 

I ^ - I and  an  approximate  solution  of  liquation  3-25  is  found  by  first 

obtaining  C."  from 
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I 


, "I  Jk:i<  - ' V1  1 «^2/(^a)| 


and  then  ns  ini; 


where  i n i t i a 1 1 v 


T //  2 "R 

L,t:  ’ - a 


,k‘=R  + lc 


1/2 


(3-28 ) 


(3-2'.)) 


le  ll  + l2 


(3-30) 


R ^ ( 2 n - 1 ) 


Both  t and  R are  varied  in  an  iterative  solution.  t is  modified 
«.  n e 

to  account  for  the  second-order  term  in  liquation  2-20,  and  R is  modi- 

n 

t i ed  to  account  tor  variation  in  1 ^ with  decreasing  conduct i v i ty . 

In  the  iteration  procedure,  and  R arc  represented  in  subsequent 
steps  in  the  iteration  by 

si/-1 


*c  = ‘l  + J'l2  + ■’  *3 


(f 


R 


I 5-32 ) 


I2n  - I)  - j .|)  I o 1 

«r  p p 


(3-53) 


where  ’ and  C are  values  from  the  previous  iteration  and  I ) 

I P ..r  p p 

is  evaluated  using  liquation  3-2"  with  values  of  t and  q determined 

by  c"  and  a , respect  ivelv. 

P P 

Case  lb.  Re(-t)  > 1,  horizontal  polarization. 

The  solution  for  horizontal  polarization  is 


'■K  - > 


1 + IziV/lalO 

R n 


,k“-R  + V 


( 5-5  1 I 


where  R = 2 n and  is  relatively  independent  of  conduct  ivitv , and  i 
is  determined  as  described  above. 
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Case  2a.  | t | 1,  vertical  polarization. 

lor  waveguide  upper  boundaries  determined  by  normal  or  weakly 
disturbed  nighttime  boundaries,  a graz  i ng-earth  incident  angle  approx  i 
mat  ion  is  often  appl  icable  lor  lower  mode  numbers.  It  the  magnitude 


of  t is  small  and 

the  conductivity  is  greater  than  about  10  mho/rn. 

1 ^ 

j gr 

explj  /(> 

I exp  [-.i  1 ( - t )'V/  “j 

( 3-35  1 

Making  the  furt be 

r approximation  that 

■> 

t:“  • - 

a 

(3-36) 

which  should  be  true  it  t 

is  small,  the  solution  to  1 

quation  3-25  can 

be  written  as 

l 2 

K (:Ll 

n ' .*> 

it?)' : 

1/2  , 1 

) ' s?) 

4 = 

s> 

1 

' ■ ■ 3 2 

- t- 

1/2 

1 

' ti  1 

(3-42) 

H i s defined  bv 
> n 

1 

. . 
n 

IJ11  - 5 

i -n  1 

gr 

(3-38) 

1 1 

l n i t 1 a 1 1 v , R is 

n 

tinen t ly  modif  ied 

defined  f 

in  an  itc 

or  high-conductivity  ground  and  is  subsc- 
■ration  loop  to  account  for  low  conductivity 

Cast  2b.  t 1,  horizontal  polarization, 
lor  these  conditions, 


T 4 3/2l 

I J exp  (j5 ■'/<>!  exp  |^-j  -1-4*  J 


(3-39) 
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I lie  solution  for  H is  the  same  as  for  vertical  polarization  except 


that  K is  replaced  lu- 
ll 1 


R = ( 12n  - I )/(, 


and  does  not  vary  with  conductivity. 


(3-10) 


base  .>a.  Relt)  1,  vertical  polarization. 

"hen  the  ionospheric  reflection  altitude  is  h i «h , there  result 
what  is  known  as  a "detached  whispering  gallery"  mode.  Here 


„r  J exp  ( j / 2)  expj^-j  ' (-t  )‘^2J 


(3  — 4 1 | 


C~  = X“  < 1 - \“ 


v is  equa  1 to  the  value  of  C“  that  is  obtained  from  I 


" ‘y*  . u | 

■(?)  a ,'f . sa 


with  R replaced  hv 
n ' 


("V  ') 


quation  3-37  but 


( 3-43  ) 


(la so  3b . Re  ( t I 1 , 

I he  same  equations  used  for  vertical  polarization  are  used 
for  this  case. 

I II.  RATION  ALGORITHMS . liquation  3-15  is  solved  using  Newton's 
iteration  method.  Given  an  initial  guess,  t , an  improved  guess  is 
obtained  from 


t , t + At 
nl  n() 


I 3 Hi 


where 
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I 

I 


I 


I 


11 


\l  t 

„)B|t  ) - 1 

nU  n() 

'"'no 

'A,,' 

-ilqT  - t + y I B ( t 
i n()  o nO 

[M'nO 

I v\ ' ( t - > ) + q.K  ,(t  - V 

[ J n()  <>  i J nO  o 

i 3 - 1 . 


\n  iteration  procedure  using  Equations  3-44  and  3-45  and  an  initial 
ness  from  an  appropriate  approximation  to  liquation  3-25  is  used  to 
make  t small,  for  Newton's  method  to  converge  to  the  appropriate 
'line  for  a specified  mode,  a sufficiently  pood  guess  is  required  for 

1 ill)  ’ 

\ second  iteration  procedure  is  now  used  to  obtain  the  exact 
eigenangles  of  liquation  3-1.  The  improved  guess  determined  above  is 
used  with  the  Newton-Raphson  and  anisotropic  ionospheric  reflection 
coefficient  algorithms.  Ihe  procedure  is  to  compute  a correction  term 


n 


-l;(") 


Ml  | • 


0=0 


( 3-  In 


where  is  the  value  for  the  present  iteration.  A new  value  is  used 
o 

defined  using 


and  Iquation  3- In  is  applied  iteratively  until  AO  becomes  less  than 
some  preset  criterion. 

VLF  FIELD  STRENGTH  EQUATIONS 

The  111  DCOM  code  provides  two  alternative  procedures  for  com- 
puting the  received  electric  field  strength  produced  by  an  arbitrarily 
oriented  transmitting  dipole.  Ihe  first  method,  the  Kkll  approximation, 
assumes  that  the  ionospheric  profile  and  ground  conductivity  vary 
slowly  as  the  signal  propagates  through  the  earth- ionosphere  waveguide. 
As  long  as  there  are  no  abrupt  transitions  and  the  total  transition  is 
not  severe,  this  first-order  approximation  is  adequate  and  an  easily 
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interpretable  set  of  computations  is  provided.  On  the  other  hand,  the 
IV KB  approximation  becomes  of  questionable  use  when  sharp  changes  in 
the  ionosphere  or  ground  characteristics  occur.  Implicit  to  the  WKB 
approximation  is  that  each  mode  can  be  traced  as  the  signal  propagates 
along  the  waveguide.  Note  from  the  example  of  figure  3-3  the  difficulty 
of  mode  identification  if  there  is  a sharp  transition  from  a daytime 
to  a nighttime  profile.  In  addition,  a sharp  transition  would  scatter 
energy  from  one  incident  mode  into  other  modes,  violating  another  VV K B 
assumption  that  all  the  energy  would  remain  in  the  same  mode. 

These  problems  do  not  occur  with  the  second  method,  the  Mode 
Conversion  procedure.  \t  a transition  region,  mode  identification 
need  not  occur  for  this  procedure,  and  energy  scattering  into  other 
modes  is  considered,  fxcept  for  special  cases,  however,  the  mode  con- 
version computations  are  not  as  easily  i nterprctablc  as  is  the  IVKB 
met  hod . 

Both  procedures  segment  the  waveguide  into  N cascaded,  uniform 
waveguides  where  the  propagation  characteristics  do  not  change  in  any 
segment.  Thus,  the  eigenangle  computations  defined  earlier  arc  appli- 
cable for  either  technique.  Before  procoding  to  a detailed  description 
of  these  procedures,  the  hi. DCOM  transmitted  power  input  options  will 
be  discussed.  The  input  quantity  is  the  effective  frce-space  antenna 
input  power,  lor  a lossless  short-wire  antenna  in  fre't'  space,  the 
input  power  and  current  moment  are  related  by 


I 


efs 


l<f(,fl 


kb 


13-  48) 


where 

f = frequency  (kill) 

I = current  moment  (ampere  meters) 
m 

However,  it  is  the  equivalent  orthogonal  (vertical,  broadside,  and 
end fire)  dipole  moments  that  arc  used  in  IVHPCOM: 
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»»■ 


P - P . cos 

V CIS 


> ■> 

P,  - P . . si  n*“  • cos**  { 
I c t s 

■)  ■> 

I’,  = I’  ..  sin4".'  sin"! 
H els 


( 3-49 l 


where  (m.t  I igure  1 - 1 J 

■ = dipole  zenith  angle  I deg  I 

angle  between  the  .\-a.\is  and  the  dipole  a.xis  fdegj  . 

WKB  Approximation 

lhe  KkB  formulation  used  in  Isi  I'fOM  is  basically  that  described 
in  Reference  .3-5  and  is  repeated  here  for  completeness.  The  three 
orthogonal  dipole  components  defined  above  are  considered  a vertical 
component  and  two  horizontal  components,  one  aligned  with  the  direc- 
tion of  propagat  ion  iciidfire  and  one  at  right  angles  to  the  direction 
of  propagation  (broadside  . three  components  ot  the  electric  field 
are  computed:  I , I , and  I . the  basic  equations  used  for  a nonuni- 
form ionosphere  for  each  mode  are  given  below: 

\ 

I - 1 g ’’  V/m  , ( -5-50.) 

xn  n ■ n • n 

where  n is  the  mode  number,  and 


l x , y , 


( 3 - 5 1 j 


nKzn 


"I/-’ 


B-  it  vn 


■ g • (5-52) 
I . ■ n xn  ( 


11.054  1 1 

T) 


exp  [ - j kd ( S - 1| 


(5-551 


f frequency  I kHz) 

a earth's  radius  (km) 

d path  length  I km  I 

k wavenumber  (km-') 

(\klt  approximation  to  the  sine  of  the  cigcnanglc, 
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1 


I'he  's  are  representative  excitation  factors  computed  as  a geometric 
mean : 

, N 


ii  ill! 


I A - a ! 


I lie  superscripts  "1"  and  "N  " indicate  values  at  the  transmitter  and 
receiver,  respectively.  The  g's  are  height-gain  factors  evaluated 
at  the  respective  transmitter  and  receiver  antenna  altitudes.  Both 
's  and  g's  will  he  defined  on  the  following  pages. 

lo  account  for  nonuniformities  in  ground  and  ionospheric 
parameters  along  the  propagation  path,  two  approximations  are  made. 

An  average  value  is  computed  for  the  factor  S in  liquation  5-53. 

There  are  N computation  points  for  a trunk,  evenly  spaced  at  distance 
d|,  so  that 


■in  • j + sin  6^ 


N -1 
P 


S|"  I 


^ + V*  sin  0. 


i = 2 


I JK 

I J 


(5-55) 


where 


. is  the  eigenanglc  at  the  ith  computation  point. 

Ihc  average  attenuation  rate  and  relative  phase  velocity  are 
computed  from 


i = -8.(>9k  I m(S)  x 1 o'"*  dB/Mm 
. = RelSJ  , 


(5-5(i ) 


(3-57 | 


wlie  re 


c = velocity  of  light  (hm/s) 
v = phase  velocity  (km/s)  . 

Mode  Conversion  Method 

The  Wl.DCOM  mode  conversion  algorithm  is  an  adaptation  of  the 
.V)SC  model  (Reference  5-4J.  This  model  assumes  that  there  is  wave- 
guide invariance  normal  to  the  waveguide,  and  that  the  reflections 
from  the  slab  transitions  may  he  neglected.  This  is  consistent  with 
Wait  (Reference  5 - 5 J , who  provides  an  argument  demonstrating  that 
neglect  of  reflections  does  not  introduce  much  error.  However,  this 
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aspect  should  he  explored  further  for  very  sharp  transitional  regions 

found  in  some  nuclear  environments. 

Die  mode  conversion  algorithm  formulation  initially  assumes 

that  a unit  amplitude  wave  exists  in  each  mode  in  t lie  transmitter 

slab.  lhe  quantities  \l'  relate  the  energy  from  an  incident  mode  k 

K t li 

in  the  transmitter  slab  to  the  j mode  in  slab  p.  The  mode  conver- 
sion coefficients  are  computed  by  enforcing  the  boundary  condition 
at  each  slab  interface.)  that  all  tangential  field  components  must 
be  continuous. 

Hie  three  orthogonal  electric  field  components  I (■  = x,y,  and 
are  computed.  Mien  the  waveguide  is  uniform,  the  component  compu- 
tation- are  similar  to  the  1\KB  solution;  ic,  for  mode  n 


. 1 

n ■ n ■ n 


(3-581 


I i tin.  "to  .eiK-ral  nonuniform  waveguide  case, 


y </  p 

/ 111,11 


g 1 \ 1 , 

■ n ' in  mn 


(3-59) 


0 . 0514 1 


r f/a  r 

sin (d/a) 


f- j k J sinlU^J  + 


N d -| 

s,,,,  lj  ;>  - 


(.5-60) 


d(,  individual  slab  length  (km)  . (3-61 

Np 

lhe  mode  conversion  coet 1 ic icnt , A , includes  the  eltects  ot  mode 

mil 

conversion  and  propagat ion  in  all  slabs  prior  to  the  receiver  slab. 

Np 

Derivation  ot  A ' is  performed  in  Reference  3-6.  lhe  Kf.IK.6M 
mn 

version  is  a simple  extension  of  this  derivation  to  consider  cases 
where  the  number  of  eigcnangles  is  not  necessarily  the  same  for  all 
slabs.  The  mode  conversion  coefficients  are  determined  by  first 

■J 

solving  for  the  parameter  a“  from 
* 1 m , n 


47 


t'\ 

O L 

- ua  cos  p yj  i - 

(d-(i'.l) 

t'.’  - 
1 

t1'.  - kh/k 
o i a 

(3-70J 

1 he  IV  s 

a rc 

A i rv 

functions  and  the 

1 ' s will  be  de f i ned 

in  the 

nex  t 

subsec t 

i on . 

1 he 

a^  are  related 

m , n 

to  the  mode  conversion  coeff 

’ ic i ent s 

In 

in , n 

s i n ( ''  ) 
p m 

‘ m , n , 2 

s l n ( " ) 
n 

( 3 - ~ 1 J 

the  at. 
j K 

and 

aP+1 
l ,n 

(right  sides  of  1 

.quat  ions  a -(id  and  .5 — ( > 

5 | a re 

so  1 ved 

using  f 

rout 

's  L-V 

decomposition  al 

gori thm  for  the  solut 

ion  of 

s i mu  1 - 

taneous  linear  equations  ( Reference  3-7J  of  the  form 


| A | X - B . (3-7JJ 


11k  matrix  |.\J  is  or  order  ( p- 1 ) by  (p-1),  and  \ and  B are  ( p - 1 
ordered  vectors.  As  |A|  and  B are  known,  the  vector  X can  be  solved 
for.  An  improved  estimate  for  X is  obtained  by  first  defining  an 
error  residue, 

K [ A | X - B [A|6X  . (3-"3) 


A correction  to  X is  made  for  6X : 
o 

X = X - vx  . 


( o-7  l ) 


o 

Ibis  procedure  is  performed  iteratively  until  'X  becomes  negligible. 

The  above  procedure  for  computing  A*  is  not  used  where  the 
propagation  characteristics  between  two  slabs  are  almost  identical, 
lor  this  case, 


A1’ 
in , n 


m ,n 


- ik 


c* 


s i n ( * 


,p-i 

m 


H 


( 3-75} 
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VLF  Anistropic  Excitation  Factors 

Hu‘  A' s in  liquation  3-52  are  excitation  Factors,  where  the 
t i rst  subscript  indicates  the  dipole  component  and  the  second  subscript 
indicates  the  electric  Field  component.  They  are  computed  From 
(References  5-5  and  5-9) 

AVX  = \-Js 

vv  ~ ^ ii  KjlM  * iiKit)(I  f iUi ) 


= B — 


I,  = -Av:/S‘ 


1 1 * II  ^ II  1 ( 1 ~ l 1 l^l  1 


'l  l = J n'V  I + „R„  II  1 * 1R1) 


- -Av=/S 


A Bx  = ABz/S 


B H + iM"'1  - 


M * n n ^ m ^ 


in- 


B= 


where 


- siK„(l  + iRi  I ( 1 + MkJ  , 


5/  1 

- S' ' “/at-  ( " j /do  | o=o 


= eigenangle  (sin  O^J 


(5-76 J 


F(.o)  = mode  equation  (see  liquation  5-lj  . 

I he  It's  are  ionospheric  reflection  coefficients  described  in 
Section  2,  and  It's  are  ground  reflection  coefficients  defined  in 
Appendix  C. 


50 


4 


VLF  Anisotropic  Height  Gain  Factors 

the  g's  in  F(|iiation  3 - 52  arc  height  gain  factors  (Reference 
3 Si  and  are  computed  from 


g_/a  - c 


) 


l\-,(t  )F.  - IV,  ft  )!•'  J k/k 


"x  it 


J 


‘Wi  “ Ki(V|;: 


I . - Kj  (t .)!■ ' , 

T2  VVVV  1 i 


1 3-78  | 


where 


g_  = e 


(“) 


It)  I - lv  ( t ) F 


‘Wl  ' lVi  (thJI  2 


'i  W - 'W 

I , w'(t  J - QW  ( t I 
J o Jo 

'.3  = ‘VV  - VW 

1 I = h2<V  - VVV 


(3-79) 

(3-80  I 
(3-81  ) 
(3-82) 
(3-83) 


~ is  the  transmitter  or  receiver  altitude  (km),  ~ is  the  bottom 

hot 

of  the  ionosphere  (ic,  the  altitude  where  the  reflection  coefficient 
integration  is  stopped),  and  the  quantities  Q and  Qj j are  defined  in 
Append i x C . 

(he  Airy  function  arguments  t,  t and  t|  arc  defined  by 


t = k“(C~  - 2(1)  - z ) / a J 


t = k *•((:“  + — ) 
o a a 


t, 


k”  (C* 

a 


hot 


J/aj 


(3-84  j 
(3-85) 
(3—80 ) 


where  (.'  is  the  cosine  of  the  eigenangle,  and  I)  is  the  altitude  where 
the  modified  index  of  refraction  is  unity.  lor  the  mode  conversion 
computations  I)  and  z.  are  assumed  to  he  zero. 
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VLF  Isotropic  Height  Gain  and 
Excitation  Factors 


Uhcn  magnetic  field  etlect  . are  unimportant,  the  isotropic 
ionosphere  approximation  is  used  and  the  field  ec|uations  are  simpli- 
fied. From  I qua t i on  3-76  when  1R||  and  nU  are  zero, 

vy  = 'Bz  = lv  = 0 (a-8"| 

I hr  exc  it.it  ion  and  height  gain  factors  required  to  compute 

tile  and  \ components  of  the  electric  field  are  X X \ 

v:’  I :’  A By  ’ 'x  ’ 

f-v  * a,Ki  3 • In  the  Airy  function  notation  of  IVa  i t (Reference  3-d), 
thev  are 


•a.  a '■:>  l 


t - >•„  - ■ 


i 


1 1 1 - c|i\ , 1 1 1 


\ . 

1 V 

1 

i 

w ' ( t 

: - 

-V  + 

i|i 

l\-,(t  - 

>vr 

( 3-88  ) 

where 

the  q 

’ s and  i| 

■ ' s a re  era 

l 

1 uat 

ed  for 

vert  i 

cal 

polar 

■ i z 

at i on , 

and 

\ is  a 

cons 

taut.  \ 

is  era  lu 
By 

at  ed 

from 

the  si 

line 

expres 

si  on  but 

with 

the  q’ 

s and 

Vs  cv 

aluated  for 

hor 

izontal  pola 

ir  i z 

at  l on  . 

1 n fqli 

at  i on 

3-88  a 

11  ionospheri 

c -dependent 

quant i tie 

s are 

computed 

at 

the  t 

rans- 

• i tt er 

and 

rcce i ver 

1 oca  t i on  s . 

To 

be  applied 

d i r 

ect  J v 

i n 

the  e 

1 cc  t r i c 

f i el d 

eqliat 

ions  (Equations  3-5 

2 o r 

5-55 

and  3- 

56  J 

, t he 

exc i tat  ion 

f ic  t or 

defined  by  1, 

quat  ion  3-S.s  requires 

that 

I he  height  gain  factors  are  computed  from 

<’.Sj  J Wj  (t  - z)  l\'(t  ql\,it  ij 

- h'  i ( t - y J • l\ ' ( t ) i|l\ ^ ( t ij  , 

where  \ i the  normalized  transmitter  or  receiver  altitude: 


(3-90 i 


1 

( 


I 

I 


! 


lor  vertical  polarization  g^  0 and,  likewise,  for  horizontal  polari- 
zation g (l.  Also,  for  vertical  polari  at  ion. 


■ h 1 1 - y . ) + B( t ] IV j ( t - y __ ) 

Kz  i~  U , It  - yj  + Bit  J l\  ( It  - y . ) 


SIGNAL-TO-NOISE  RATIO  CALCULATION 

Lor  the  vertical  electric  tield,  the  \ 1.1  sijjnal-to-noi.se  ratio 

in  dB  is 


s\U  dll  log...  (I  /.\  ) 
Jo  z o 


(3-yd  i 


I .rs 


jy/. 


10 


I l'ani/-0  I 


I 5-95  : 


where 

frequency  (khz) 

1 Bandwidth  (llzj 

noise  figure  idB  above  275  • Boltzmann's  constant), 
i defined  i;  Reference  5-10,  which  also  contains  values  for  a 

a;. 

rounded  vertical  dipole.  The  values  are  not  valid  for  other  polari - 


.at i on  . 


ELI  EIGENANGLE 

Only  one  mode  propagates  at  1.1.1  , and  it  lias  been  found  that 
almost  any  reasonable  initial  estimate  for  tlu  eigenangle  will  result 
in  convii . i nee . lu  compute  the  initial  value,  a typical  eigenangle, 

( , is  first  defined.  \ first -order  improvement  is  made  to  this 
value  a 1 Glow--.  I he  vertical  ground  reflection  coefficient  is 
assumed  unit)  and  the  modi-  solut ion  is  approximated  h\ 

(.„!<„  - 1J  0 . (3-9-1) 

further,  tin  assumption  is  made  that 
(t>  i • 

,,  R , e ■ , 1 3-95  j 

where  ■ i a complex  constant.  The  reflection  coefficient  is  computed 
u in.;.  a typical  incident  angle  . from  liquations  5-1)  I and  T-95  it 
follows  that  the  initial  value,  .,  can  he  found  from 


‘f 
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4 


0. 

1 


I 

I 


(3-96) 


| 

I 

I 


-1 

= COS 


I 

I 


I tt  COS  !) 


I he  exact  value  of  the  III  eigenangle  is  determined  by  the  Newton- 
Raphson  iteration  technique  defined  by  Equations  5-46  and  5-17  using 
as  an  initial  value. 

l 

ELF  FIELD  STRENGTH  EQUATIONS 

tor  transmitted  field  strength  at  l.l. I'  frequencies,  the  verti- 
cal l-tield  radiated  from  a surface  horizontal  dipole  is  computed. 
Since  there  is  only  one  mode,  the  field  is  computed  from 


( 3 - H7  I 


where  superscripts  I and  K refer  to  the  transmitter  and  receiver  loca- 
tions respectively  and  the  excitation  factors  are  defined  as  they  are 
tor  \ 1 .1  cases,  and  the  height  gain  factors  are  computed  from 


C e-*kC: 


A * 


ikCz 


i 


( 3 - 1)8  J 


1 + „K 
a a 


i Mi z t - i kCz 
c + i l'i  • 


8,. 


(3-99) 


1 ♦ ,11, 


g_ 


.ikCz  . - ikCz 

e + K e 

M II 


1 ♦ , K 

ii  ii 


(5- l no i 


these  are  flat  -earth  approximations  which  are  valid  at  11.1  where  the 
ionospheric  reflection  coefficients  are  integrated  to  ground  (ie, 

2bot  = • 

Ihe  hklt  procedure  (described  earlier  for  VI. I field  strength 
computations)  is  used  for  the  I l.l  cases.  The  ionospheric  ionization 
profile  is  analyzed  at  two  additional  locations  (edges  of  Iresnel  zone 
at  midpath)  to  determine  if  the  ionosphere  is  sufficiently  uniform  to 
assume  the  correctness  of  the  WKB  approximation.  The  MTKlOM  output  is 
flagged  it  if  is  determined  that  the  IV K It  approximation  might  be  invalid 
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SECTION  4 


SKYWAVE  AND  GkUUNDWAVE  CALCULATIONS  FOR  LF 


INTRODUCTION 

In  the  II  program,  computations  are  performed  for  the  I._  (ver- 
tical], f ^ (perpendicular),  and  L (end)  field  components.  Mien  an 
arbitrarily  oriented  dipole  antenna  is  used,  all  three  orthogonal 
components  from  each  of  the  three  orthogonal  dipole  moments  are 
summed . 


the  II  propagation  model  uses  ray  theory  with  correction  terms 
for  diffraction  around  the  surface  of  the  earth.  The  ray  geometry  is 
adjusted  for  a nonuniform  ionospheric  disturbance.  Elevated  antennas 
arc  modeled,  hut  a built-in  model  simplification  (used  to  increase 
computation  efficiency)  is  not  applicable  for  antenna  altitudes  higher 
than  aircraft  altitudes  (ie,  b-8  km). 

LF  FIELD  STRENGTH  EQUATIONS 


Hie  total  field  strength  for  signals  propagated  at  IT  is  rep- 
resented as  the  sum  of  a ground wave  and  a number  of  skywaves.  Speci- 
t’ical  1>  , 


M 

L 1 

m-  1 


,m 


14-1  i 


T = I. 


M 

IX 

111=  1 * 


,111 


(4-2) 


where 


m 

» s 

m=  1 


I , m , 

x 


l 4-3) 


I - total  vertical  field  strength  (V/m) 
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! 


i 


l 


I total  field  strength  perpend  i cular  to  direction  of 

propagation  (V/m) 

I total  field  strength  in  direction  of  propagation  (V/m) 

1 vertical  component  of  groundwave  (V/m) 

1 horizontal  component  of  groundwave  (V/m) 

I vertical  field  strength  for  skvwave  with  m ionospheric 

z , m 

reflections  (V/m) 

I perpendicular  field  strength  for  a skvwave  with  m iono- 

y , m 1 ' 

spheric  reflections  (V/m) 

I end  field  strength  for  a skvwave  with  m ionospheric 

x , m • ' 

reflections  (V/m). 

When  the  ionosphere  is  anisotropic,  the  skywaves  for  vertical 
and  horizontal  polarization  are  coupled.  This  point  will  he  discussed 
in  detail  later. 

Radiated  Power 

the  fields  are  formulated  in  terms  of  radiated  power.  In  the 
computer  model  the  effective  frec-spacc  lossless  dipole  antenna  input 
power,  I’  . , is  specified.  I he  three  orthogonal  (vertical,  broadside, 
and  end  fire,  i dipoh  powers  are  used  in  the  1.1  model  computations: 


'\ 

P COS 

C’ t s 

■) 

(4-4) 

PB  = 

1 ' e t s S‘n"’ 

• cos"; 

( 4 - 5 J 

I',  I’  . sin“  cos";  , 
i els 

where  (see  figure  1-1) 

■ = dipole  zenith  angle  (deg 

; = angle  between  the  X-axis  and  dipole  axis  (degj . 


(4-t>) 


Groundwave  Field 

The  groundwave  field  is  given  by  (References  4-1  and  4-2) 

o.si-yv 

I.  = : — V/m  (vertical  dipole)  (4-7) 

zg  d 


o.3Py2w 


V/m  (horizontal  dipole 

cndf i re ) 


V/m  , (horizontal  dipole 

broadside) 


0 = conductivity  at  transmitter  (mhos/m) 

j - relative  dielectric  constant  at  transmitter 
= permittivity  of  tree  space  (farads/mj 
= wave  frequency  (radians/s) 

j = sT- 1 , 

and  where  d is  the  transmitter-receiver  separation  distance  in  kilo- 


I lie  loss  functions  W and  IV  arc  defined  by  (Reference  4-2) 


o-jXts  lVl(ts--V  "l  (ts  - 


1 TS  - ^ 


W,  ft  ) 


h . ( t i 
1 s 


(1-101 


and  Wv  is  expressed  the  same  way,  except  that  Q_  is  replaced  by  Q in 


I .qua  t i on  1-10  and 


«■  (*)  V). 

e ' 


(4-11) 


where 


k = wavenumber  (km  ) 
a = equivalent  earth's  radius  (km) 
and  t (related  to  the  eigenangle)  is  defined  below. 

In  the  LI-  calculations,  a^  is  defined  as  4/a  the  actual  earth 
radius  to  correct  for  tropospheric  refraction.  The  Q's  are  defined 


Q,  = 


- 1 


C U! 

o 


1/2 


o + i > ■ a 

k ' g o 


(4-12) 


Q = 


(4-13) 


where  and  arc  average  conduct  ivitv  and  relative  dielectric 

g g 

constant  values  over  the  entire  path. 

W j is  the  \ir>  function  of  the  fir  t kind  (Reference  4-3J. 
I he  arguments  y , y,  are  defined  b\ 


(4-1 4 j 


(4-15) 


where  h and  h are  transmitter  and  receiver  antenna  heights  in  kilo- 
t r re 

meters.  The  values  of  t are  roots  of  the  differential  equation 
(Reference  4-1) 


dt 

-dlf  = 


(4-lb) 


and  q is  either  (for  lq_)  or  Q (for  i;  ) . 

The  values  of  t.  are  found  by  integrating  the  differential 
liquation  4-lb  with  respect  to  q using  a fourth-order  Runge-Kutta 
integration  technique.  Solutions  for  the  starting  condition  (q=0) 
are  known. 

When  Q (Q^  or  is  small  (less  than  one),  the  integration 

is  performed  by  integrating  hquation  4-lb  from  0 to  Q using  as  a 
starting  value 


t ( q = 0 ) 
s 


- - jn/3 


(4-17) 


‘ 
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where  Z is  the  sth  zero  ot  \ . . The  quantity  A.(-Z)  is  the  eon- 
s i i 


ventionally  defined  \irv  function  (Reference  1-4).  The  first  10  zeros 
are  shown  in  Table  4-1.  When  more  than  10  terms  are  required  in  the 
integration,  higher  order  zeros  arc  computed  from  (Reference  4-4) 


( y . ) 


2/3 


35 


181223 


4Sv' 


288  v 
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207360v 


1 S(>8357  1 ' 
2441()0v8  . 


1 lid 


.■> 

8 


(4s  - 3) 


(4-18) 
(4-19  I 


After  defining  the  initial  values  of  t.,  the  integration  is 
performed  using 

' (P,  + 21’  + 2f>-  + I’,  | , (4-20) 


t (q.  . ) = t (q.  ) 
S 1 + 1 s 1 | 


Table  4-1.  Zeros  of  A.(-,)  and  A^(-a). 


Sequence 

Number 

ex?  = Zero  of 

Ai  (-.*) 

u'c  = Zero  of 

A; (-a) 

1 

-2.33810741 

-1.01879297 

2 

-4.08794944 

-3.24819758 

3 

-5.52055983 

-4.82009921 

4 

-6.78670809 

-6.16330736 

5 

-7.94413359 

-7.37217726 

6 

-9.02265805 

-8.48848673 

7 

-10.04617434 

-9.53544905 

8 

-1 1.00852430 

-10.52766040 

9 

-11.93601556 

-11.47505663 

10 

-12.82877675 

-12.38478837 

where 

Aq  = qj+,  - q.  (1-21) 

= * I M : » ts(q  )JAq  (1-22) 

I’  = f[q.  + 1/2, '.q,  t.(q  ) + l/2l’.|,\q  (4-23) 


(4-24 j 
(4-25J 


>\3  = 

•y. 

< 

i 

'4 

t'lHj  + Aq,  ts(q.J  + 

Ihe  solution 

is  started  using 

l 

"'■'j  - 

v 1/21’ I Ah 
1 Aq  • 

(4-26) 


Solutions  for  t arc  found  using  Aquations  l - 1 1>  through  4-26 
followed  by  defining  the  sth  term  in  the  summation  in  liquation  l Jo. 
Terms  are  det  i tied  tint  il  the  last  term  is  less  than  10  \ 

for  vertieal  polarisation  and  very  low  ground  conductivity 
or  for  horizontal  polarisation  in  all  cases,  Q lie  Q_  or  Q | becomes 
large.  I he  procedure  for  determining  t is  similar  to  that  used  for 
small  values  of  t).  Ilu  substitution 

'-k 


(4-27) 


is  made  in  Aquation  l-l(i,  resulting  in 


dt 

dp 
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_J 

- v 
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I lie  integration  is  performed  by  taking  steps  in  p instead  of 
q in  Aquation  l - It),  and  the  starting  values  for  p=(>  are  now  defined  by 

,,  - jn/3 

ts(p=0)  = a.  e J , (.)-2ii) 

where  the  u ' s are  the  zeros  of  A. |-u'.  I he  first  Id  zeros  are  riven 

s i 

in  'I able  4-1.  When  more  terms  are  needed,  higher  order  zeros  are 
defined  by  (Reference  4-4) 
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Skywave  Field  Strength 

I'lie  total  skywave  field  strength  is  computed  as  the  sum  of 
field  strengths  of  up  to  five  skywaves.  Figure  1-1  shows  the  geometric 
idealisations  used  in  performing  the  skywave  calculations.  The  rays 
bounce  hack  and  forth  between  the  earth  and  an  idealized  sharply 
bounded  ionosphere.  The  ionospheric  boundary  is  defined  at  the  alti- 
tude, hj, , that  is  in  the  center  of  a region  in  the  ionosphere  where 
most  energy  is  reflected  (see  liquation  1-1  lb  for  the  definition  of 

V- 

Path  1 in  Figure  1-1  shows  a case  where  the  receiver  is  beyond 
the  geometric  horizon  for  a one-hop  path.  In  this,  instance  the  ray 
is  assumed  to  graze  the  earth's  surface,  and  energy  propagates  the 
remaining  distance  by  diffraction  around  the  curved  earth.  Ibis 
situation  is  referred  to  as  propagation  beyond  the  caustic.  In  the 
idealization  used  here,  the  diffracted  distance,  Dj , is  assumed  to  be 
uniformly  divided  between  transmitter  and  receiver  ends. 

Paths  2 and  3 illustrate  ray  paths  that  arc  within  the  geome- 
tric limit  and  show  the  geometry  for  a nonuniform  ionosphere.  The 
angle  and  distance  quantities  identified  on  the  figure  are  used  later 
in  the  discussion  of  the  iteration  procedures  that  define  the  geometry. 

Fields  for  skywaves  with  N,  N + 1,  and  N + 2 ionospheric  re- 
flections are  always  calculated,  where  N is  the  minimum  number  of  sky- 
wave  hops  required  to  traverse  the  path  without  diffract  around  the 
earth’s  curvature  are  also  considered.  If  N is  3 or  greater,  skywaves 
with  N - 1 and  N - 2 hops  arc  considered.  If  N is  2,  the  field  for 
a skywave  with  one  hop  is  computed. 

The  general  equations  for  describing  the  skywave  field 
strengths  are  given  first.  Then  the  procedures  defining  the  skywave 
geometry  and  the  various  parameters  that  are  used  in  the  general 
equations  are  provided. 

CFNFRAI.  I 1 1 i 1 .1 ) liQlJAT  I ONS . Because  the  model  is  developed  to 
perform  calculations  when  the  ionosphere  is  nonuniform  and  anisotropic, 
the  equations  defining  the  multihop  field  strengths  arc  cumbersome. 
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1 he  following  general  nomenclature  is  used  to  simplify  writing  the 


equat i ons : 


Ionospheric  reflection  coefficient 
for  the  ith  hop  of  a j-hop  path 


Antenna  foreground  factor  at  the 
transmitter  for  a j-hop  path 


Antenna  foreground  factor  at  the 
receiver  for  a i -hop  path 


l’o lari  tat  i on 


vertical  to  vertical 
horizontal  to  vertical 

vertical  to  horizontal 

horizontal  to  horizontal 

vertical  field, 
vertical  dipole- 

horizontal  field, 
broadside  dipole 

vertical  field, 
endfi re  d i pole 

vertical  field, 
vertical  dipole 

horizontal  field, 
broadside  dipole 

vert ical  field , 
end  fire  dipole 


Ground  reflection  coefficient  of 
diffractive  correction  factor  for 
the  ith  hop  of  a i-hop  path 


vert ical 


hor i zonta 1 


Vertical  field  for  a j-hop  shywave 
Horizontal  field  for  a j-hop  shywave 
Convergence  coefficient  for  a j-hop  shywave 


Inverse  distance  field  constants  that  account  for  the  radiated 


power  and  total  path  distance  arc  first  computed.  These  are 


o. isrj'- 

I = 1 I vert  i ca  1 ) , or 

zo  d 


(end  fire)  (4-32) 


0.  1.31' 


where  P , l’j  , 1’^,  and  d are  as  defined  previously. 

Che  field  for  a one-hop  skywave  is  computed  by  first  coni[)uting 
the  intermediate  quant i t y 


„W,V\  , JIV.HB 

1 1 zoRl 1 I 1 + voRll  T1 


.,1111.  1IB  . .All  AX 
>vl  voRl  1 11  + ::okl  1 11 


(.4-34  ; 


(4-35) 


where  \ \ for  a vertical  dipole  or  X 1 for  an  endtire  horizontal 

d i po  1 e . 

Chen  the  field  is  computed  from 


' zi  [zi,'iaai 
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(4-57) 


The  field  for  a two-hop  skywave  is  computed  by  first  computing 


the  intermediate  quantity: 


(4-38 J 
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The  field  for  ;i  greater  number  of  hops  is  computed  by  a continuation 
of  the  process  described  above. 

A vector  sum  and  mis  field  sum  arc  computed.  The  vector  sum 
is  obtained  from 


11  = 5*2 

L 1 

n=L 


V/m 


(4-14  j 


(n  .\+J  \ 

•:’«  • & 'if  "*  • < i-.s, 

In  liquations  4-44  and  4-45,  L is  the  maximum  of  1 and  N - 2. 

The  rms  sum  is  computed  because  of  uncertainty  in  the  absolute 
phase  of  each  skywave  in  a disturbed  ionosphere,  and  because  the  pre- 
diction of  deep  nulls  or  high  peaks  in  field  strength  calculations 
for  a single  path  distance  may  he  misleading.  Ivi-IXIOM  outputs  the  rms 
quantities  when  the  field  strength  is  to  be  computed  only  at  the 
receiver.  The  vector  sum  is  output  when  the  field  strength  is  computed 
at  intermediate  points  along  the  path. 

1 . 1.l . VATl-.l)  ANihNNA  CONS  IDf.KAT  IONS.  When  elevated  antennas  are 
used,  up  to  four  paths  are  considered  for  each  number  of  hops  as  illus- 
trate... in  figure  4-2.  The  field  strengths  resulting  from  propagation 
via  paths  a^  , a,,  a^,  and  b,  are  added  vectorial ly  to  obtain  a compo- 
site field  strength  for  the  particular  skywave.  The  special  definition 
of  the  antenna  foreground  factors  for  the  case  with  elevated  antennas 
is  discussed  later  in  this  section.  The  four  paths  shown  in  figure 
4-2  are  sufficient  to  define  propagation  for  elevated  antennas  at 
less  than  (>-7  km.  for  higher  altitudes,  four  more  individual  paths 
need  to  be  considered. 


I NTHITMl.l)  I ATT!  CALCULATIONS.  The  skywave  geometry,  ionospheric 
reflection  coefficients,  convergence  coefficients,  foreground  factors, 
and  ground  reflection  coc ‘‘i c i cut s or  diffraction  loss  factors  are 
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Figure  4-2.  Geometry  for  four  paths  for 

elevated  transmitters  or  receivers. 

computed  by  performing  the  following  sequential  steps: 

1.  Search  procedure  for  defining  the  skywave  geometry  and 
locating  the  ionospheric  and  ground  reflection  points. 

2.  Interpolation  procedure  for  defining  ionospheric  and 
ground  electrical  properties  at  the  reflection  points. 

is.  Computation  of  ionospheric  reflection  coefficients. 

I.  Computation  of  convergence  coefficients. 

5.  Computation  of  ground  reflection  coefficients  or  terms 
defining  diffractive  losses  around  the  earth's  bulge 
and  antenna  foreground  factors. 


facli  of'  these  steps  is  described  below: 


Search  for  Reflection  Points.  Re  fleet  ion  coefficients  and 
re  fleet  ion  points  (altitude  and  position  along  path)  are  determined 
for  skywaves  with  m hops.  Sky ’waves  with  in  less  than  N (where  \ is 
the  minimum  number  of  hops  needed  to  reach  the  receiver  with  no 
diffraction,  are  assumed  to  diffract  over  part  of  the  path.  An  exam- 
ple i:  shown  in  figure  4-1,  where  the  skywave  with  one  reflection 
(path  li  is  assumed  to  propagate  with  one  maximum-distance  hop  (graz- 
ing incidence  at  the  earth)  over  part  of  the  distance  between  trans- 
mitter and  receiver  and  to  propagate  the  remaining  distance  U’j)  hv 
diffraction.  In  the  example  shown,  two  hops  arc  the  minimum  number 
of  hops  that  will  reach  from  transmitter  to  receiver  with  no  diffrac- 
tion. (Iround-based  transmitters  and  receivers  are  shown  in  the  example 
for  simplicity.  Modifications  to  account  for  elevated  antennas  are 
given  later. 

1 he  steps  in  locating  the  reflection  points  are  as  follows: 

1.  A reference  altitude  and  an  ionospheric  scale  height  are 
determined  at  specified  intervals  between  transmitter  and 
receiver.  The  reference  altitude  and  ionospheric  scale 
height  are  described  below  under  interpolation  Procedure. 

Skywaves  with  less  than  \ hops  are  assumed  to  have  graz- 
in'.’ incidence  ( ; = PO  degrees).  I he  distance  between 

g 

transmitter  and  receiver  is  divided  into  m equal  intervals, 
and  the  reflection  point  of  each  hop  is  assumed  to  occur 
at  the  center  of  each  interval.  A trial  reflection  alti- 
tude for  each  hop,  h , is  picked  (a  starting  value  lower 
than  the  expected  reflection  altitude],  and  the  angle  of 
incidence  with  the  ionosphere,  {• . , is  computed.  By 

interpolation  between  values  calculated  at  the  specified 

intervals,  the  actual  reflection  altitude,  h , is  found. 

a 

If  h^  is  greater  than  h^  , then  h is  increased  and  the 
process  is  repeated  until  h^  is  approximately  equal  to 
h . 
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;S.  lor  the  r:i > path  with  \ hops,  a trial  angle  ( i • calcu- 
lated, a nil  ionosphere  reflect  ion  point*  and  ground  reflec- 
tion points  are  found.  Ihe  aiu.lt  i i then  varied  until 
tin  Nt  h ground  reflection  point  coincide-  w 1 1 ti  the  receiv- 
er. In  carrying  out  this  process,  it  is  assumed  that  all 
ground  reflection  angles  are  equa 1 to  the  take  off  angle; 
however,  the  reflection  altitude  and  hop  distance  max  In- 
different for  each  hop. 

1.  Stop  5 is  repeated  for  N + 1 and  \ + 2 hops  to  define 
the  remaining  ray  paths.  Ihe  dashed  line  for  path  A in 
I igure  1-1  shows  one  step  in  the  interactive  process  in 
which  i is  varied  to  find  the  location  of  the  reflection 
points  with  three  hops  that  reach  the  receiver. 

Interpolat ion  Procedure.  Ionospheric  parameters  are  evaluated 
at  the  trial  reflection  point  hy  an  interpolation  procedure.  The 
variables  used  in  the  interpolat ion  procedure  are  shown  in  I igure  1-d. 
Ionospheric  parameters  are  known  at  distances  dj  and  d,  that  bracket 
tin  trial  distance  d, . A linear  interpolation  procedure  is  used  to 
define  electron  and  positive  ion  vertical  profiles  at  the  distance 
d . 

t 

Ihe  ionospheric  parameters  used  in  the  interpolation  procedure 
are  a reference  altitude  and  an  ionospheric  scale  height.  Both  are 
defined  in  terms  of  BlhJ,  the  imaginary  part  of  the  ionospheric  index 
of  reflection.  I he  ionospheric  scale  height,  h^_,  is  defined  by 

ii  B{h  J/dT,  [H(hJ]  ’ < 

I hi  reference  altitude  hR  is  defined  by  determining  the  altitude  for 
wh  i c h 

B(hRJ  = 0.0*1  . (4-47 ) 

I hi  reference  altitude  h„  and  scale  heights  h are  obtained  at 

K s 

distance  d by  linear  interpolation  between  the  values  at  dj  and  dv 
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Figure  4-3.  Variables 

The  actual  reflection  altitude, 

B (h  i = V~2  cos"  d> . 

r 


in  interpolation  procedure. 

h , is  then  estimated  from 
r 


( 4-48  l 


where  <{>.  is  the  ionospheric  ang  1 e of  incidence  shown  in  figure  4-3. 
Using  i uuation  1-45  and  the  definition  of  h..  and  h results  in 

K S 
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I' 


V" 


yf ~ '■’Os“  <p. 


o.  04 


(4-4‘J) 


lhe  ionospheric  region  centered  about  the  altitude  h defined 
by  Iquation  4-41)  has  been  shown  to  be  the  region  from  which  most  V 1.1 
and  LI-  energy  is  reflected.  Strictly  speaking,  it  is  applicable  for 
an  isotropic  ionosphere  and  therefore  for  daytime  or  disturbed  condi- 
tions. However,  because  of  the  steep  vertical  gradient  of  the  iono- 
spheric profile  for  normal  nighttime  conditions,  the  reflection 


altitude  is  not  sensitive  to  the  criterion  used,  and  therefore  liqua- 
tion 1-19  is  used  for  all  conditions. 

When  elevated  antennas  are  used,  four  ray  paths  are  considered 


as  shown 

in  figure  4- 

4. 

In 

the  iteration  procedure, 

the  great -circle 

d i stance 

cove  red 

by  c 

ach 

ray 

path  is 
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dbl 

= "hi 
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+ Rbl 
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= V 

> 

*1)2 

* Ub2  • 

(4-55 j 

In  all  cases  ( x is  either  a or  b), 

I'  = a 1/2-  v . - f ■ ) , 
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(4-56 1 

where 


-./  Ue  \ 

'xi  = S,n  la  +h~r°SYxi 

. -1  / “e  \ . . 

‘xi  \ae  + htr  / XJ 

. -1  / Ue  \ 

u . = s,n  mrr,s  v 

XJ  \ e rc  / 


I*.,.  = sin 

1 x i 


■ u . ~ sin 


(4-57) 


(4-58) 


(1-59) 


and  a is  the  equivalent  earth's  radius. 

Ionospheric  Reflect  ion  Coefficient.  I he  above  search  proce- 
dure results  in  the  definition  of  the  ionospheric  incident  angle  <}k 
and  the  distance  from  the  reflection  point  to  the  transmitter,  d,. 
Linear  interpolation  is  used  to  define  electron  and  ion  density 
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profiles,  magnetic  propagation  azimuth,  dip  angle,  and  the  magnetic 
field  strength  at  the  reflection  point.  With  these  quantities  deter- 
mined, the  ionospheric  reflection  coefficients  are  obtained  using 
procedures  described  in  Section 

For  elevated  antennas,  reflection  coefficients  are  computed 

for  onlv  two  incident  angles,  1 , and  , (see  figure  4-1).  Reflec- 

al  bl 

t ion  coefficients  for  path  a ^ are  used  for  path  a,,  and  reflection 
coefficients  for  path  bj  are  used  for  path  b,. 

V onv  ergence  factors . When  the  distance  between  transmitter 
and  receiver  is  less  than  the  distance  to  the  caustic  and  a uniform, 
concentric  ionosphere  is  assumed,  the  convergence  factor  for  an  m-hop 
ray  is  (Reference  1-5) 


J j ' 

/ a + h , ,,  1 11  + ;l  - a l-'os  C 4 /2m J 1 

I _c a l Jm  sin  I ■ / -in  i a e e 

\ a / sin  . J I ( h ( + a ) cos  I, /2ml  - a 


( 4-()(l  J 


ionospheric  height  (km) 
number  of  hops 
d / a^,  ( rad i ans  ) 

ground  range  between  receiver  and  transmitter  (km)  . 


for  a nonuni torrn  ionosphere  but  neglecting  any  tilt  of  the 
ionosphere  (ie,  the  height  varies  but  tiie  ionosphere  is  locally  con- 
centric!, it  can  be  shown  (see  Reference  1<>;  that 


( 1-(«1  J 


where 


bill 


the  reflection  altitude  at  the  1th  reflection  point  (km 


i 


whore  d is  the  distance  to  the  caustic  for  m hops, 
cm 

liquation  4-63  is  applicable  for  a uniform  ionosphere.  for  a 
nonuniform  ionosphere,  liquation  1-63  is  used  with  the  approximation 
that 

m 

h = 1 Y\  h (1 ) , (4-67 j 

a m a 

i = 1 

where  h ( ( l J is  the  reflection  altitude  for  the  1th  hop . 

The  convergence  factors  for  path  a^  are  used  for  path  av  and 
coefficients  for  path  b ^ are  used  for  path  b?. 

1 1 round  deflect  ion  tlocff ic i cut , foreground  factors,  and  Diffrac- 
tion Loss  terms.  The  influence  of  the  ground  on  the  propagated  signal 
depends  on  ray  path  geometry.  The  geometric  idealizations  that  are 
used  in  computing  the  influence  of  the  ground  at  the  transmitter  loca- 
tion are  defined  in  figure  4-5. 


Figure  4-5.  Geometry  to  compute  the  influence  of 
the  ground  at  the  transmitter. 

When  the  geometric  optics  limit  applies,  the  foreground 
and  a,  are 

+ R e -'*'■’’1  (vertical  dipole)  (4-68) 

RV 


factors  for  paths  a. 

„VV 


1 i 


= 1 
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( 1-83  ) 


S,  = Sh, 


V 


who  re 


Sbj  tot. il  r;i>  path  length  tor  path  1> j 

Sit,  total  ray  path  length  (or  patli  b,  . 

Similar  definitions  are  used  tor  the  receiver  foreground 

, , - i K S i . - iRS  > 

factors,  except  that  the  phase  correction  terms  e • and  e • “ are 

not  used  since  the  total  phase  correction  is  included  in  the  trans- 
in i 1 1 e r t e rm . 

When  the  antenna  altitude  is  very  small  or  when  the  elevation 
angle  at  the  ground  is  near  zero,  paths  a^,  a,,  b^  , and  h,  merge  and 
the  foreground  factors  for  the  composite  path  are  computed  from 


U 


lilt 
' T.i 

, V I 


»1>V 


K (v  v.) 
v ' I 1 


" r. 1 y,l  ’ -v  i J 


(vertical  dipole)  (1-84 J 


(horizontal  dipole-  (4-85) 
broads i de) 

(horizontal  dipole-  (4-86) 
endf i re) 


g wt 

lht  loss  functions  W and  W are  evaluated  using  liquation  1-1()  with 

vl/3 


1 


V . 

• 1 


K) 

($f 


(4-87) 


(4-88 J 


who  re 


li  = transmitter  altitude 
t r 

h.  ionosphere  reflection  altitude  for  the  ret  lection 
point  nearest  the  transmitter. 

1 he  switch  from  using  liquations  1-68  through  4-83  to  using 

s 1-81  through  1-88  is  made  when  i j 0.03  radians. 

figure  1-6  shows  the  geometry  used  tor  calculating  ref  lection 


losses  or  diffraction  losses  at  intermediate  ground  reflection  points. 
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c)  g;ov-  o»'.  >•  ■•on  <b>  kycnd  GrCMriinc  optics  «gion 

Figure  4-6.  Geometry  to  compute  the  influence  of  the 

ground  at  an  intermediate  ground  reflection  point. 

Figure  4-6(a)  illustrates  the  case  where  the  receiver  is  well  inside 
the  geometric  limit.  Here  the  losses  arc  defined  by  the  Fresnel 
reflection  coefficients  given  in  liquations  4-72  and  4-73.  lor  the 
situation  depicted  in  Figure  4-6(1)),  the  loss  is  computed  by  evalu- 
ating liquation  4-8,  with  Vj  and  y defined  using  h.j  and  h.,,  and  X 
defined  using  liquation  4-11  and  by  the  distance  between  adjacent 
ionospheric  reflection  points. 
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APPENDIX  A 

DEFINITION  OF  THE  S MATRIX 


The  elements  of  the  partitioned  S matrix  (liquations  2 -3a 
through  2 -3d ) are  defined  in  terms  of  an  intermediate  T matrix.  The 
defining  equations  are  [Reference  2-5) 


I J + 'll 


14/C:  * ‘"'ll 


,/C  - CT 


12/t:  + ‘42 

12/t:  - ‘ 4 2 

31  + T34/C 
-,/c  - T-, 


+ «32/(; 


- r32/C 


(i  = cos 


5 sin 


= complex  angle  of  incidence  (deg) 

and  the  intermediate  matrix  I relates  the  electromagnetic  field  com- 
ponents to  their  derivatives.  In  turn,  T is  defined  in  terms  of  the 
susceptibility  matrix  M: 
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I in;ill>,  the  suscept  i hi  1 i t > M matrix  is  defined  ;i: 
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S I II 


collision  frequency  (s  ) 

wave  frequency  (s  ) 

free  space  permcab i 1 i t v Ihenrvs/m) 
o 

■ magnetic  dip  angle  Idcg) 

a earth's  radius  (km) 

1)  altitude  where  index  of  refraction  is  unity  (km) 

4 direction  of  propagation  east  ot  magnet  ic  north  (degj 
- altitude  (km)  . 

The  matrix  M described  above  is  general  in  that  it  represents 
a single  ionospheric  constituent.  three  species  are  considered  by 
this  model:  electrons  and  positive  and  negative  ions,  thus  for  each 
species  there  will  be  a matrix  M , and  the  matrix  elements  needed  for 
computing  the  1 matrix  elements  arc  found  from 

M , fM.  - 

i = l 


APPENDIX  B 


STARTING  SOLUTION  FOR  THE 
ANISOTROPIC  INTEGRATION  TECHNIQUE 


Closed-form  expressions*  for  starting  values  for  the  reflec- 
tion coefficient  integration  are  obtained  by  first  finding  the  ti%o 
upgoing  have  solutions,  q and  q , , of  the  Booker  quart ic 

q ' + lb_q^  + (>b_q2  + lb | q +b()  0 . (B-l) 

The  parameter  q is  equal  to  the  product  of  the  index  of  refraction 
and  tin.  cosine  of  the  angle  of  refraction  in  the  ionosphere,  and  the 
coefficients  are  functions  of  the  susceptibility  matrix  M: 

b_  = B _ / ( 4 B . J 

a a 1 

i),  = K,/  (<>B  j J 
I'l  = BJ  / (4  B j ) 

bn  = VB4 

B4  = 11  + Ma.V> 

B_  - S(M  + M ) 

a 1 a a 1 

B,  - (C“  + M__)  ( 1 + M . ) + M M_ . 

J a a II  I a a I 

-<i  * *!,-){(:“  * M„J  * M,5M-, 
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- M..(C2  + M , , ) M - . - (1  + M..INL.M.,  - M.-M^lC  + M._) 

la  22  a 1 11  2a  a2  12  21  aa 


'From  Reference  2-a. 
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I' he  K matrix  elements  are  then  determined  from 
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APPENDIX  C 

GROUND  REFLECTION  COEFFICIENTS 


The  symbols  R represent  the  ground  reflection  coefficients 
referenced  to  the  altitude  where  the  ionospheric  reflection  coeffi- 
cient computations  are  terminated.  They  are  defined  as  follows: 
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In  Equations  C-3  through  C-(>, 

(.'  ; cosine  of  the  eigenangle 


n“  = 1 - 2 |D 


hot 


j/a  , 


when 


I)  - reference  altitude  where  the  index  of  refraction 
is  unity 

= altitude  where  ionospheric  reflection  coefficient 
calculations  are  terminated. 


hot 


The  double  subscripted  A's  are  defined  in  terms  of  Airy 
functions,  W.  and  W,,  of  the  I irst  and  second  kind  as  defined  by  Wait 
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in  Reference  3-2.  In  l.quat ions  C-7  through  C-14,  primes  indicate 
derivatives  with  respect  to  the  argument. 
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